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Phonons are essential for understanding the thermal properties in monolayer tran¬ 
sition metal dichalcogenides, which limit their thermal performance for potential 
applications. We investigate the lattice dynamics and thermodynamic properties of 
M0S2, MoSe2, and WS2 by first principles calculations. The obtained phonon fre¬ 
quencies and thermal conductivities agree well with the measurements. Our results 
show that the thermal conductivity of M0S2 is highest among the three materials 
due to its much lower average atomic mass. We also discuss the competition between 
mass effect, interatomic bonding and anharmonic vibrations in determining the ther¬ 
mal conductivity of WS2. Strong covalent W-S bonding and low anharmonicity in 
WS2 are found to be crucial in understanding its much higher thermal conductivity 
compared to MoSe2. 


I. INTRODUCTION 

Monolayer transition metal dichalcogenides, MX2 (M = Mo, W; X = S, Se), have aroused 
much interest recently due to their remarkable properties for applications in next-generation 
nanoelectronic and energy conversion devices^^-. Since all these applications are closely 
related to its thermal properties, it is necessary to investigate the lattice dynamics and ther¬ 
modynamic properties of MX2. For instance, high-performance electronic devices strongly 
depend on high thermal conductivity for highly efficient heat dissipation, while low thermal 
conductivity is preferred in thermoelectric application. 

In semiconductors, heat is carried by the atomic vibrations that are quantized as 
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phonons^. Theoretical predictions based on the phonon Boltzmann transport eqnation have 
fonnd that monolayer WS 2 has the highest thermal condnctivity (142 W/mK) at room 
temperatnre, then followed by M 0 S 2 (103 W/mK) and MoSe 2 (54 W/mK}^. However, the 
measured thermal conductivity for monolayer M 0 S 2 and WS 2 is 34.5 ± 4 W/mK- and 32 
W/mK-, respectively, which is much lower than the theoretical predictions. Furthermore, 
various phonon properties such as interatomic bonding and anharmonic vibrations still lack 
a unihed understanding. The parameters that affect phonon transport include crystal struc¬ 
ture, atomic mass, interatomic bonding, and anharmonicit}*^'— . Generally there are four 
rules for hnding a nonmetallic crystal with higher thermal conductivity, including (i) lower 
average atomic mass, (ii) stronger interatomic bonding, (iii) simpler crystal structure, and 
(iv) lower anharmonicity^. All monolayer MX 2 compounds have similar crystal structures, 
while conditions (i) and (ii) imply a larger Debye temperature, and condition (iv) means 
smaller Griineisen parameter. Recent theoretical investigation has provided a quantitative 
analysis of the roles of mass, structure, and bond strength in thermal expansion and ther¬ 
momechanics of MX^^. However, the roles of mass, interatomic bonding, and anharmonic 
vibrations in phonon transport still remain uninvestigated. Glear knowledge of the underly¬ 
ing physics will be helpful for understanding and modulating the thermal transport in MX 2 , 
for example, through doping other M or X atoms^^ii^. 

Here we investigate fundamental vibrational properties to understand thermal trans¬ 
port in M 0 S 2 , MoSe 2 , and WS 2 . The measured phonon frequencies are well reproduced in 
our calculations. The thermodynamic properties are calculated within quasi-harmonic ap¬ 
proximation, and the calculated thermal conductivities agree well with the measurements. 
Gombining first principles calculations and the Slack model, the roles of mass, interatomic 
bonding, and anharmonicity in thermal transport are clearly revealed. 


II. METHODOLOGY 

All calculations are implemented in the Vienna ab initio simulation package (VASP) 
based on the density functional theory (DFT) method^. The Perdew, Burke, and Ernzerhof 
(PBE) parametrization within the generalized gradient approximation (GGA) is used for the 
exchange-correlation functional. A plane-wave basis set is employed with the kinetic energy 
cutoff of 600 eV. A 15x15x1 k-mesh is used during structural relaxation for the unit cell 
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TABLE I: Calculated lattice parameters and band gap of monolayer MX 2 .Experimental data are 
also given in parentheses for comparison. 


M 0 S 2 

MoSe2 

WS 2 

a (A) 3.165 (3.160 

3.300 (3.288 

3.163 (3.154 '^) 

Eg (eV) 1.81 (1.88 A 

1.56 (1.57 

1.97 (1.95 A 


^ Referencel^ 

Reference!^ 

Reference's 

ReferenceSi 

Reference's 

until the energy differences are converged within 10“® eV, with a Hellman-Feynman force 
convergence threshold of 10“^ eV/A. We maintain the interlayer vacuum spacing larger than 
12 A to eliminate the interaction with periodic boundary condition. 

In the calculation of phonon dispersion, the harmonic interatomic force constants (IFCs) 
are obtained by density functional perturbation theory (DFPT) using the supercell approach, 
which calculates the dynamical matrix through the linear response of electron density^. 
The 5x5x1 supercell with 5x5x1 k-mesh is used to ensure the convergence. The phonon 
dispersion is obtained using the Phonopy code with the harmonic IFCs as input^S. 


III. RESULTS AND DISCUSSION 

A. Crystal structures 

Monolayer MX 2 has honeycomb structure with space group P6m2i^ as shown in £g. [H An 
M atom layer is sandwiched between two X atom layers, connected by covalent bonds. The 
optimized lattice parameters of all studied MX 2 are shown in table [H Our GGA calculations 
overestimate the lattice parameters by 0.16%, 0.36%, and 0.29%, respectively, which is a 
general feature of the GGA functional. 

The electronic structures of all studied MX 2 are calculated by DFT method. As shown 
in table [H the calculated band gap is consistent with the measurements^^—. The total and 
atom projected density of states (DOS) are shown in £g. El The total DOS from -7 to 10 
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eV is mainly composed of M-d and X-p states, and the bands on each side of the band gap 
originate primarily from M-d states, which is in agreement with previous work—. Due to a 
less localized DOS of W atoms, the overlap between W-d and S-p state in the valence band 
of WS 2 from -7 to 0 eV is larger than other two materials , indicating a strong covalent p — d 
bonding. Similar large overlap between W-d and S-p state in WS 2 from 0 to 12 eV tends to 
increase the widths of the conduction bands. 

Fig. |3] presents the electronic charge density of all studied MX 2 in the [110] plane. For 
M 0 S 2 and MoSe 2 , the highest charge density is found to be on the Mo atoms due to the 
strongly localized DOS of Mo atoms, while for WS 2 , S atoms have the highest charge density. 
As shown in £g. El the W-S bonding is the strongest, whereas the Mo-Se bonding is weaker 
than the Mo-S bonding, which is consistent with the projected DOS in £g. [21 Usually, strong 
interatomic bonding and low average atomic mass imply a large Debye temperature, leading 
to a high thermal conductivity. From £g. |2l it is found that the d state of transition metal 
M play an important role in determining the interatomic bonding, which will further affect 
the heat transport in these three materials. 


B. Phonon spectra 


The phonon spectra of all studied MX 2 structures are calculated using the supercell ap¬ 
proach, with the real-space force-constants calculated in the density-functional perturbation 
theory (DFPT)^ within the Phonopy code^^. Fig. 0] presents the phonon spectra along sev¬ 
eral high symmetry directions, together with the corresponding projected phonon density 
of states (PDOS). The primitive cell of monolayer MX 2 contains 3 atoms, corresponding to 
three acoustic and six optical phonon branches. The average acoustic Debye temperature 
for monolayer MX 2 is determined frontal 


1 

where dj for each acoustic branch i [i 
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where H is Planck constant. 
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0Ji,max is the phonon frequency at the zone boundary of the 


Ath acoustic mode, and ks is Boltzmann constant. The calculated Debye temperature On 
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for M 0 S 2 , MoSe 2 , and WS 2 are 262.3 K, 177.6 K, and 213.6 K, respectively, which is in 
good agreement with previons resnlts, i.e. 260-320 K for M 0 S 2 estimated from specihc-heat 
measnrement^, 197.3±6.6 K for MoSe 2 estimated from photocondnctivity measnrements^, 
210 K for WS 2 estimated from the Lindemann formnlaj^^. 

Concerning thermal vibrations and the bonding forces, the Debye temperatnre is a mea- 
snre of the temperatnre above which all modes begin to be excited and below which modes 
begin to be frozen ont^^. We hrst investigate the relation between vibrational modes and 
the Debye temperatnre by calcnlating the projected PDOSs for the M(XY), M(Z), X(XY), 
and X(Z) vibrations in MX 2 stnctnres as shown in £g. 01 Similar to the diatomic linear 
chain model, the scale of the aconstic (optical) phonon branch is dominated by atoms with 
larger (smaller) mass in three materials. As the mass ratio of all stndied MX 2 (mM/'m 2 x) in 
table mi show, the aconstic phonon vibration in the PDOS is governed by the larger mass. 
The mass ratio of M 0 S 2 is most close to 1, while that of WS 2 is mnch larger 1. Therefore 
the low-freqnency acoustic phonon branches of M 0 S 2 up to 233.9 cm“^ are mainly from 
the Mo(XY), Mo(Z) and S(XY) vibrations due to similar mass, whereas those of WS 2 up 
to 182.3 cm“^ are mainly from the W(XY) and W(Z) vibrations due to the much larger 
mass of W atoms. In contrast to other two materials, the mass of transition metal atoms 
in MoSe 2 is smaller than the mass of chalcogenide atoms. Thus, although all Mo(XY), 
Mo(Z), Se(XY), and Se(Z) vibrations contribute signihcantly to the low-frequency branches 
of MoSe 2 up to 157.5 cm“^, the PDOS of the Se(XY) and Se(Z) vibrations is higher than 
other two vibrations due to the relatively larger mass of Se atoms. 

Furthermore, low average atomic mass M, besides the strong interatomic bonding men¬ 
tioned above, can lead to a large Debye temperature as well. The mass differences in M 0 S 2 , 
MoSe 2 , and WS 2 are also shown in table m For M 0 S 2 and MoSe 2 which have similar bond¬ 
ing characteristics as shown in fig. [3l the M of M 0 S 2 is approximately two thirds the M 
of MoSe 2 , and the calculated 6d is about 1.5 times larger. For WS 2 , although it has close 
mass to MoSe 2 , the strong covalent W-S bonding as shown in £g. [3l can result in relatively 
larger Debye temperature as implied in previous work-. Therefore, although the average 
atomic mass plays a key role in determining the Debye temperature, the effect of inter¬ 
atomic bonding can not be neglected and should be taken into account when studying the 
thermal transport properties of monolayer MX 2 . 

In addition, we investigate the vibrational mode of all studied MX 2 . Since the monolayer 
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TABLE II: The mass ratio average atomic mass M,aiid Debye temperature 6d of all 

studied MX 2 . 


Structure mM/nr 2 X M (amu) 

Od (K) 

M 0 S 2 

1.50 

53.36 

262.3 

MoSe2 

0.61 

84.63 

177.6 

WS 2 

2.87 

82.66 

213.6 


TABLE III: Theoretical determined optical phonon frequencies (cm at the T point. The exper¬ 
imental results are also given in parentheses for comparison. 


Structure E” 



E' 


^'1 


4' 


M 0 S 2 

278.4 

(283 

a) 

377.2 

(385 

) 397.7 

(404 

) 461.6 

(470 '^) 

MoSe2 

163.3 

(167 


235.7 

(240 

) 279.3 

(282 ^ 

) 345.2 

(351 ") 

WS 2 

291.6 

(298 

a) 

350.8 

(357 

) 413.0 

(419 

) 432.8 

(438 

^ Reference^ 


Reference^ 

References^ 

References^ 

MX 2 belongs to the D^h point group, the optical lattice-vibration modes at F point can be 
thus decomposed as 


Topucai = dl"(IR) + A;(R) + E'(IR + R) + E"(R), (3) 

where IR and R denote infrared- and Raman-active mode respectively. Table HTTI lists the op¬ 
tical phonon frequencies at the F point. The calculated phonon frequencies are in agreement 
with the experimental results, and the discrepancy is less than 3%. The LO/TO splitting is 
very small and can be neglected hereSSi^S. 

The schematic vibrations for the phonon modes are shown in £g. [5], where one A'^ and 
one E' are acoustic modes, the other A 2 {E') are IR (both IR and R) active as shown in 
Eq. ([3]). A 2 and A'^ modes vibrate along the 2 ;-direction, and E' and E” modes vibrate in 
the x — y direction. As shown before in fig. 01 in the case of E'(LA/TA) in monolayer M 0 S 2 , 
the Mo and S atoms vibrate with similar amplitudes; for A 2 (ZA) in M 0 S 2 , the vibrations 
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of Mo atoms have much larger amplitudes. For both i7'(LA/TA) and A 2 (ZA) in monolayer 
MoSe 2 , the Se atoms vibrate with greater amplitudes than Mo atoms. The vibration of W 
atoms dominates S atoms in both i7'(LA/TA) and A 2 (ZA) vibrational modes in monolayer 
WS 2 due to the large mw/'m 2 S- 

C. Thermodynamic properties 

The Griineisen parameter 7 , which describes the thermal expansion of a crystal on its 
vibrational properties, provides information on the anharmonic interactions. The larger the 
Griineisen parameter indicates the stronger anharmonic vibrations. The expression for the 
Griineisen parameter is given by^>^ 


7 = 


3aBVm 


a 


V 


(4) 


where a is the linear thermal expansion coefficient, B is the bulk modulus, Vm is the molar 
volume, and Cy is the isometric heat capacity. 

Table HVl compares the calculated isometric heat capacity, bulk modulus, and linear ther¬ 
mal expansion coefficient with experimental results at 300 K. The isometric heat capacity 
can be calculated as 


f9E\ 1 g^n(q)/fcBr 

\~k^} (e'-4q)ABT_i)2’ 


(5) 


where T is temperature, and a;n(q) is the phonon frequency of the n-th branch with wave 
vector q. The calculated values of Cy for M 0 S 2 , MoSe 2 , and WS 2 at room temperature are in 
good agreement with the experimental results^^— . The bulk modulus B and linear thermal 
expansion coefficient a are calculated using the quasi-harmonic approximation (QHA), which 
takes the hrst-order anharmonicity into account^. The obtained B for M 0 S 2 is in agreement 
with the experimental valued. For MoSe 2 , the computed B is fallen in the range of measured 
values^*^. The great discrepancies between the calculated and experimental B are seen for 
WS 2 . Since the bulk modulus is used to describe the stiffness of MX^^^, the bonding in 
monolayer WS 2 is found to be much stiffer comparing to other two materials. The calculated 
a for M 0 S 2 and MoSe 2 at 300 K is is fallen in the range of measured values^^'^'^, while 
larger than the measured a for WS 2 . 






TABLE IV: Comparison between the calculated and measured isometric heat capacity Cy (J 
mol“^ bulk modulus B (GPa), and linear thermal expansion coefficient a (10“® K“^). The 

experimental results are also given in parentheses. 


Structure Cy 

B 

a 


M 0 S 2 

62.97 (63.55 

52.3 (53.4±1.0 ’’) 

17.4 (10.7 

(82 <^) 

MoSe2 

68.75 (68.60 

57.3 (45.7±0.3 0 (62 §) 

19.5 (7.24 

(105 

WS 2 

63.49 (63.82±0.32 

') 77.9 (61±1 j) 

14.8 (6.35 



^ References^ 

ReferenceSS 
Reference's 
Reference's 
® References! 

^ Reference's 
s Reference^! 

^ Reference's 
^ ReferenceSS 
j References! 

^ Reference's 

As shown in fig. | 6 l the temperatnre-dependent Griineisen parameter is calculated using 
Eq. dl]). The Griineisen parameter can also be calculated by averaging the mode Griineisen 
parameter 7 n(q), 


mode 

lave 


1 


X]7n(q)Gy,„(q), 

n,q 


( 6 ) 


where Gy^n(q) is the mode heat capacity. The mode Griineisen parameter is given by 


7n(q) = 


ao dun{q) 


(7) 


(:n„(q) da 

where Oq is the equilibrium lattice constant at 0 K. Fig. EJalso shows the calculated 7 ^°'^®, 
which is consistent with the Griineisen parameter calculated using Eq. (|1]). 

The frequency dependence of mode Griineisen parameter of M 0 S 2 in the irreducible BZ 
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is plotted in fig. [TJ and the mode Griineisen parameter along symmetry directions is shown 
in the inset. Similar to diamond, graphite and graphene^, negative 7 za is observed at 
low frequencies (under 50 cm“^ in M 0 S 2 ). The average Griineisen parameter is negative at 
low temperatures because only the low-frequency phonons are excited, as shown in £g. | 6 l 
Our results are in agreement with previous calculations^*^. The calculated at room 

temperature is 1 . 22 , 1 . 20 , and 1.15 for M 0 S 2 , MoSe 2 , and WS 2 , respectively, indicating that 
M 0 S 2 has the strongest bonding anharmonicity among three materials, while WS 2 has the 
weakest anharmonic vibrations. The bonding anharmonicity obtained from analysing the 7 
of all studied MX 2 is consistent with previous work—. 

According to Slack’s expression^*^^, assuming that only the acoustic phonon modes par¬ 
ticipate in the heat conduction process, the lattice thermal conductivity in the temperature 
range where three-phonon scattering is dominant, can be given as 


°2r ^ 


( 8 ) 


where M is the average mass of an atom in the crystal, 5^ is the volume per atom, n is the 
number of atoms in the primitive unit cell, and A is a constant which is given by^ 


2.43 X 10"® 

1 - 0.514/7 + 0.228/72 


(9) 


when the units of k;, M, and 6 are taken as W/mK, amu, and A, respectively. The obtained 
lattice thermal conductivity for monolayer M 0 S 2 , MoSe 2 , and WS 2 at room temperature is 
33.6 W/mK, 17.6 W/mK, and 31.8 W/mK, respectively, which is in good agreement with 
the experimental value of 34.5±4 W/mK for monolayer MoS 2 ^, and 32 W/mK for monolayer 
WS^-^. Although there is no experimental value for monolayer MoSe 2 , our calculated ki is a 
reasonable prediction. 

The Slack’s expression attempts to normalize the effect of mass density, crystal structure, 
interatomic bonding, and anharmonicity^AS. M 0 S 2 , MoSe 2 , and WS 2 have similar crystal 
structure. The factor M6^6 in Eq. (| 8 ]) is maximized for light mass, strong bonded crystals, 
because low average atomic mass and strong interatomic bonding lead to a high 9^, and the 
6^ term dominates the behaviour-. The Debye temperature reflects the magnitude of sond 
velocity. Higher Debye temperature results in increased phonon velocities, and increased 
acoustic-phonon frequencies as mentioned above, which suppress phonon-phonon scattering 
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by decreasing phonon populations. Therefore the high thermal conductivity of M 0 S 2 is 
related to its high Debye temperature, which is due to its much lower average atomic mass. 

Recent theoretical work has predicted that monolayer WS 2 has the highest thermal con¬ 
ductivity among all studied MX 2 at room temperature due to a large frequency gap between 
its acoustic and optic phonons^, which originates in its larger mass ratio mw/''^ 2 S- Our 
results suggest that the average atomic mass plays a key role in determining the phonon dis¬ 
persion of M 0 S 2 and WS 2 , and subsequently determines the Debye temperature from which 
the phonon velocities can be estimated. In addition, although WS 2 has similar average 
atomic mass to MoSe 2 , its strong W-S bonding leads to a higher Debye temperature. Fur¬ 
thermore, small 7 of WS 2 means low anharmonicity, which also results in a higher thermal 
conductivity. Therefore the thermal conductivity of WS 2 is determined by the competition 
between high average atomic mass, strong covalent W-S bonding and low anharmonicity. 


IV. CONCLUSION 

In summary, we investigate the lattice dynamics and thermodynamic properties of M 0 S 2 , 
MoSe 2 , and WS 2 by hrst principles calculations. The obtained phonon frequencies and 
lattice thermal conductivity agree well with experimental measurements. Our calculations 
show that the thermal conductivity of M 0 S 2 is highest among the three materials due to its 
large Debye temperature, which is attributed to the lowest average atomic mass. We also 
hnd that WS 2 has stronger covalent W-S bonding and lower anharmonicity, leading to much 
higher thermal conductivity compared to MoSe 2 . 
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FIG. 1: (a) Top view and (b) side view of crystal structure of monolayer MX 2 . The primitive 
vectors are a(a,0,0), b(a/2,-\/3o/2,0), and c(0,0,c). 
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FIG. 2: Total and atom projected DOS for (a)MoS 2 , (b)MoSe 2 , and (c)WS 2 . 
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FIG. 3: Electronic charge density of (a)MoS 2 , (b)MoSe 2 , and (c)WS 2 in the [110] plane. 
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FIG. 4: Phonon spectra and projected PDOS for (a)MoS 2 , (b)MoSe 2 , and (c)WS 2 . 
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FIG. 5: Schematic phonon vibrations with different frequencies (cm 






FIG. 6: Calculated temperature-dependent Griineisen parameter for (a)MoS 2 , (b)MoSe 2 , and 
(c)WS 2 . 
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FIG. 7: Calculated mode Griineisen parameter for (a)MoS 2 , (b)MoSe 2 and (c)WS 2 with respect 
to frequencies, and with respect to wave vectors as shown in the inset. 































